Forecasting of potential anti-inflammatory targets of some immunomodulatory plants and their constituents using in vitro, molecular docking and network pharmacology-based analysis

Most synthetic immunomodulatory medications are extremely expensive, have many disadvantages and suffer from a lot of side effects. So that, introducing immunomodulatory reagents from natural sources will have great impact on drug discovery. Therefore, this study aimed to comprehend the mechanism of the immunomodulatory activity of some natural plants via network pharmacology together with molecular docking and in vitro testing. Apigenin, luteolin, diallyl trisulfide, silibinin and allicin had the highest percentage of C-T interactions while, AKT1, CASP3, PTGS2, NOS3, TP53 and MMP9 were found to be the most enriched genes. Moreover, the most enriched pathways were pathways in cancer, fluid shear stress and atherosclerosis, relaxin signaling pathway, IL-17 signaling pathway and FoxO signaling pathway. Additionally, Curcuma longa, Allium sativum, Oleu europea, Salvia officinalis, Glycyrrhiza glabra and Silybum marianum had the highest number of P-C-T-P interactions. Furthermore, molecular docking analysis of the top hit compounds against the most enriched genes revealed that silibinin had the most stabilized interactions with AKT1, CASP3 and TP53, whereas luteolin and apigenin exhibited the most stabilized interactions with AKT1, PTGS2 and TP53. In vitro anti-inflammatory and cytotoxicity testing of the highest scoring plants exhibited equivalent outcomes to those of piroxicam.

ADME screening of the database of immunomodulatory plants. The ADME properties of the phytochemical constituents were evaluated through the use of the QikProp module, which determine some of the physiochemical characteristics that indicate a compound drug-likeness. The Lipinski rule of five served as a summary of these physiochemical characteristics. According to Lipinski's rule of five, a compound with reputed biological activity is deemed active (having satisfactory absorption and/or permeation) if it has a molecular weight less than 500 Da, fewer than 10 hydrogen-bond acceptors (Hacc), fewer than 5 hydrogen-bond donors (Hdon), 10 or fewer rotatable bonds (RBN), and has a measured log P (ClogP) less than five 28 . After filtration of database, only chemicals that obeyed at least three of the aforementioned characteristics were retained. Oral bioavailability (OB) of phytoconstituents was also estimated 29 . It shows how much of a pharmacological dosage taken orally still has an unaltered effect at the therapeutic site of action. In the database, only compounds with Identification of target proteins of the selected immunomodulatory plants constituents using network pharmacology. Due to the extravagant cost of in-laboratory screening of numerous plants along with their complicated constituents, an in-silico analysis approach was tried in this study to provide a quick, effective and high throughput technique to identify the potential immunomodulatory molecular targets of the constructed phytochemicals database. A constituent-target (C-T) network dependent on 600 plant constituents and putative targets was created utilising the screening results received from the STITCH 5.0 public database in order to uncover the mechanisms of action of the chosen plant constituents and immunomodulatory-associated target genes. STITCH is a huge database that holds a big data on chemical interactions. It connects 1.5 million genes across 373 genomes with over than 68 000 compounds, including 2200 drugs, and offers information on their interactions 30 . The function of each discovered target gene and its relationship to immunomodulation were determined using the Universal Protein Resource database (UniProt). In-depth information about more than 550 000 proteins and their functions can be found in UniProt, a resource for protein sequences and functional annotation 31 . The interactions between chemicals and genes are given "combined scores" in the STITCH 5.0 database, with stronger interactions displaying higher scores. In this investigation, only substances with interaction scores greater than 0.4 were kept (Table 1).
A total of 69 nodes (37 constituents and 32 targets) and 759 edges made up the constructed constituent-target (C-T) network ( Fig. 1), with an average of 2.435 targets for each constituent illustrating the multi-target property of the investigated compounds. The C-T network revealed that a single constituent can interact with various targets, and that targets are frequently linked to a number of different components at once. Network topological analysis was used to examine the distribution of C-T interactions on the 37 ingredients, and the results showed that apigenin had the highest percentage of C-T interactions (24%) followed by luteolin (10%), diallyl trisulfide (7.34%), silibinin (7.2%), and allicin (7.1%), respectively (Fig. 2, Table S2). Table 1 displays the discovered target genes in the network. The targeted genes AKT1 (20%), CASP3 (9%), PTGS2 (9%), NOS3 (7%), TP53 (7%) and MMP9 (7%) were found to be the most enriched ones by possessing the highest combined scores and therefore the highest percentage of interactions with the network constituents, suggesting that they may be the key nodes in that network (Fig. 3A,B). Investigation of the C-T subnetwork among the 5 top scoring constituents displayed the target genes AKT1, TP53, CASP3, NOS3 and MMP9 to be the most enriched genes (Fig. 3C). It is widely recognised that AKT1 (protein kinase B α or RAC-alpha serine/threonine-protein kinase) is a crucial signaling node in controlling the activation of adaptive immune cells 32 because it controls a variety of activities such as metabolism, cell survival, growth, proliferation, and angiogenesis 33 . Additionally, it has been shown to have substantial anti-inflammatory effects through inhibiting NF-κB-mediated transcription 34 . Meanwhile, CASP3 (Caspase-3) is implicated in the activation cascade of caspases necessary for apoptosis execution 35 . Prostaglandin-endoperoxide synthase 2 (PTGS2) is an essential enzyme for inflammatory prostaglandins biosynthesis and production 36 , while tumor suppressor protein (TP53) affects the innate immune system by secreting factors that modify macrophage activity to prevent tumour development. Additionally, it inhibits some pro-inflammatory proteins like NF-κB and STAT3, to support tissue homeostasis and prevent overreacting immune systems. Moreover, nitric oxide synthase-3 (NOS-3) is hypothesised to play anti-apoptotic roles and engaged in the suppression of genes implicated in cell proliferation. It also reduces the expression of several cytokines such as INF-γ, IL-1β, IL-6 and TNF-α, in a variety of immune cells including monocytes, eosinophils and lymphocytes. The nitrosylation of several transcription factors, including JAK/ STAT and NF-κB, mediates this action 37 . On the other hand, TP53 helps to recognise non-self-antigens, which in turn triggers anti-tumor immunity through a variety of pathways [38][39][40] . Matrix metalloproteinase 9 (Mmp9) is a secreted gelatinase that participates in the remodelling of the extracellular matrix by degrading the components of the extracellular matrix. It also controls the pro-inflammatory cytokines TNF-γ and IL-1 to control the inflammatory response 41 . A protein-protein network revealed significant gene-gene interactions (Fig. 4) indicating a high extent of synergism between them.
KEGG (Kyoto Encyclopedia of Genes and Genomes) is a database tool for deriving high-level biological system functions and applications from molecular-level data, particularly from wide-ranging molecular datasets produced by genome sequencing and other high-throughput experimental techniques [42][43][44] . To understand the signaling pathways and roles of the identified target genes, KEGG pathways functional enrichment analysis was carried out (Table 2). It can be seen from the KEGG analysis (Table 2 and Fig. 5A) that the target genes interact with 40 immunomodulatory-associated pathways, with pathways in cancer having the largest number of observed genes and the smallest false discovery rate, followed by fluid shear stress and atherosclerosis, relaxin signaling pathway, IL-17 signaling pathway and FoxO signaling pathway. To further clarify the receptor-ligand interactions associated to the disease, the constituent-pathway network (Fig. 5B) was built to connect the constituents with the signaling pathway.
The PubMed database was searched for publications linking the hit compounds to the various immunomodulatory pathways in order to validate the findings from the network pharmacology-based study (Table 3). For instance, apigenin (NF-κB and COX-2 inhibitor), was found to be valuable in lupus therapy and also for preventing inflammation in other Th17-mediated inflammatory diseases like Crohn's disease, rheumatoid arthritis, and psoriasis as well as in prevention of inflammation-based tumors overexpressing COX-2 in colon and breast as it was found to inhibit the autoantigen-presenting and stimulatory activities of the APCs (antigen presenting cells), which are required for the activation and development of lupus-associated Th1 and Th17 cells as well as B cells and causes the overactive lupus APCs, B cells, and T cells to undergo apoptosis, most likely by suppressing the expression of NF-κB-regulated anti-apoptotic molecules, particularly COX-2 and c-FLIP, which are chronically www.nature.com/scientificreports/ overexpressed by the lupus immune cells [45][46][47] . In addition, luteolin was reported to inhibit the overexpression of inflammatory mediators and catabolic factors, probably via preventing NF-κB activation. Additionally, it prevents T cell growth and their AKT/mTOR signaling. Thus, luteolin is an emerging immunosuppressant as a novel mTOR inhibitor. It also ameliorates inflammation and Th1/Th2 imbalance by regulating the TLR4/NF-κB pathway [48][49][50] . Meanwhile, diallyl trisulfide (DAT) were proved to have a potential anti-inflammatory activity by possessing an antagonistic action on TLR4 (Toll-like receptor 4/nuclear factor-κB pathway) as it inhibits the production of myeloid differentiation factor 88. DAT was reported also to increase the stabilization of IκBα by preventing its phosphorylation by IκB kinase (IKK) complex, and inhibit the NF-κB transcriptional activites 51,52 . Moreover, silibinin was declared to dramatically decrease the demyelination and inflammation signs in experimental autoimmune encephalomyelitis by upregulating the release of the anti-inflammatory Th2 cytokines and downregulating the pro-inflammatory Th1 cytokines. It also showed anti-carcinogenic and anti-inflammatory properties as a result of NF-κB transcription factor suppression. Additionally, through ERβ binding, silibinin triggers apoptosis, slows proliferation, and inhibits the production of the pro-inflammatory cytokines TNF-α and IL-17 in CD4 + T cells 53,54 . Also, allicin was proclaimed to suppress the chemokines IL-8, IP-10, and MIG, as well as the release of IL-1β from intestinal epithelial cells. This action was at least partially mediated by downregulating the mRNA levels and the involved inhibition of activation of the NF-κB pathway 55 .
Target proteins of the selected plants using combined network pharmacology. After looking at the 759 C-T interaction distributions on the selected plants, a combined plant-constituent-target-pathway network was created. The cytoscape combined score of C-T interactions was used to rate the plants. Figure 6 illustrates that Curcuma longa (29%), Allium sativum (19%), Oleu europea (17%), Salvia officinalis L (16%), Glycyrrhiza glabra (10%) and Silybum marianum (9%) had the highest number of P-C-T-P interactions which would suggest that these plants have more active substances that can act as immunomodulatory agents. In reality, there are numerous reports on the anti-inflammatory action of Curcuma longa in traditional medicine. In Indian medicinal system, Curcuma longa is used in wound healing and as an anti-inflammatory and anticancer plant 56 . Meanwhile, Allium sativum has been utilized traditionally since ancient times due to its biological properties including anti-tumour, anti-inflammatory, antioxidant, tuberculosis, arthritis, rhinitis, bronchitis, and immunomodulatory activities 57,58 . Moreover, essential oil extract of Oleu europea fruit is applied topically in Italy to alleviate rheumatism and to boost the circulation 59 . Infusion of fresh leaves of O. europaea is utilised in Tunisian traditional medicine as an antioxidant, anticancer and as a remedy for various types of inflammation 60 . In addition, Salvia officinalis L has been used to cure rheumatism and inflammation in Latin American and Asian folk medicine 61 . It also, has been utilized to treat throat and skin inflammations in European traditional medicine 62 . Numerous investigations have been carried out in recent years to identify novel biological effects of Salvia officinalis L and to document its traditional usage. Numerous pharmacological actions, such as antioxidant, anticancer, anti-nociceptive, anti-inflammatory, and antimutagenic properties, have been identified by these research 63 . While Glycyrrhiza glabra has long been utilised traditionally in Japan as an antioxidant, antiinflammatory, cancer preventative, and anti-arthritic agent 64,65 . Finally, Silybum marianum also was reported to exhibit immunomodulating and antioxidant activities by scavenging free radicals, in addition to its ability to reduce inflammation and inflammation-related pain 66 . www.nature.com/scientificreports/ Gene ontology (GO) enrichment analysis for targets. Gene ontology (GO) assigns several GO concepts to a single query sequence and classifies gene products into three distinct categories: biological process, cellular component, and molecular function 67 . By importing the selected targets into the DAVID bioinformatics tools and restricting the annotations to Homo sapiens, GO enrichment analysis was performed on the targets, revealing the most enriched pathways and GO keywords with the greatest log P value and gene counts. As demonstrated in Fig. 7A, the identified targets are linked with numerous biological processes, the most enriched ones are negative regulation of apoptotic process, positive regulation of ERK1 and ERK2 cascade, canonical Wnt signaling pathway, platelet activation and positive regulation of MAP kinase activity. The most important molecular cellular components were caveola, cytosol, protein complex, mitochondrion and nucleus. Additionally, it was shown that the most enriched molecular activities were enzyme binding, kinase activity, MAP kinase activity, protein serine/threonine kinase activity and nitric-oxide synthase activity. However, functional annotations with DAVID bioinformatics tools identified 4 BBID pathway named: MAPK_signaling_cascades, insulin-signaling, T-cell_anergy and Akt-PKB_effector_of_P13K_in_vivo,in addition to 12 BIOCARTA pathways named: insulin signaling pathway, IL-2 receptor beta chain in T cell activation, IL 2 signaling pathway, Erk1/Erk2 Mapk signaling pathway, The 4-1BB-dependent immune response, PTEN dependent cell cycle arrest and apoptosis, IL 6 signaling pathway, IL 3 signaling pathway, IL12 and Stat4 dependent signaling pathway in Th1 development, WNT signaling pathway, mTOR signaling pathway and integrin signaling pathway. Furthermore, 35 KEGG pathways including pathways in cancer, TNF signaling pathway, ErbB signaling pathway, FoxO signaling path- www.nature.com/scientificreports/ way, prolactin signaling pathway, focal adhesion and PI3K-Akt signaling pathway. In addition, 7 REACTOME pathways named ROS and RNS production in phagocytes, signaling by EGFR, negative regulation of the PI3K/ AKT network, gastrin-CREB signaling pathway via PKC and MAPK, RAF/MAP kinase cascade, VEGFR2 mediated cell proliferation and CD28 dependent PI3K/Akt signaling were identified (Fig. 7B). All of these identified pathways had P-values less than or equal to 0.01, indicating a strong correlation with inflammation.
Evaluation of network pharmacology analysis. The evaluation process including the three aspectsreliability, standardization, and rationality, was carried out according to the method described previously by Li 68 . The results were found to fulfil all the required criteria settled for network pharmacology analysis evaluation.
In vitro assays of top two target proteins; Akt1 and Caspase-3. Based on network pharmacology analysis results, additional laboratory-based in vitro screening of both Akt1 and caspase-3 inhibitory activities of the four-top scoring plants; Curcuma longa, Allium sativum, Oleu europea and Salvia officinalis L, were performed utilizing colorimetric tests. The definitive goal of these in vitro investigations is to verify the network pharmacology based analysis results through assessment of some plants-targets binding affinities. Saturosporine was employed as a positive control for Akt1 inhibitory assay since it has been showed to be a strong inhibitor of Akt1 gene 69 , where Ac-DEVD-CHO, obtained from caspase 3 colorimetric assay kit supplied by Sigma, was used as a positive control for caspase-3 gene. The results demonstrated that among the screened agents, Curcuma longa was found to have the higher inhibitory activities on both Akt1 and caspase-3 assays with IC 50 of 475 µg/ ml ± 0.031 and 500 µg/ml ± 0.053, respectively. The Akt1 inhibitory activity of other tested inhibitors was in the order of Allium sativum (IC 50 = 690 µg/ml ± 0.044) followed by Oleu europea (IC 50 = 715 µg/ml ± 0.014) and finally Salvia officinalis L (IC 50 = 740 µg/ml ± 0.063), respectively. Whereas, the caspase-3 inhibitory activity of other tested inhibitors was in the order of Salvia officinalis L (IC 50 = 590 µg/ml ± 0.045) followed by Allium sativum (IC 50 = 625 µg/ml ± 0.072) and finally Oleu europea (IC 50 = 800 µg/ml ± 0.021), respectively, (Table 4).
In vitro cytotoxicity and anti-inflammatory activity assessment. As unveiled from the prior network pharmacology study, the highest scoring plants, were retrieved to be Curcuma longa, Allium sativum, Oleu europea, Salvia officinalis L, Glycyrrhiza glabra and Zingiber officinale. Despite the fact that network pharmacology is a quick and effective method for anticipating various drug targets in complex diseases, our network analysis results need to be experimentally confirmed by examining the potential bioactivity and correlativity of those plant components in preventing the production of the leading pro-inflammatory mediators that are primarily involved with the different pathways debated above, using these monitoring indices as the basis for a comprehensive evaluation of efficacies. Therefore, to estimate the safety and efficacy of Allium sativum, Glycyrrhiza glabra, Oleu europea, Curcuma longa, Salvia officinalis L and Zingiber officinale extracts, the cell cytotoxicity 50 (IC 50 ), which is the concentration of the drug necessary to reduce the viability of cells by 50%, was established for the extract and the standard anti-inflammatory drug; piroxicam utilizing MTT test. The Allium sativum, Oleu europea, Glycyrrhiza glabra and Salvia officinalis extracts demonstrated higher IC 50 values (2324, 728, 395.9 and 127.1 μg/mL, respectively) than that of piroxicam (100 μg/ mL) illustrating that these extracts are less toxic than piroxicam (Fig. 8A). Then, using WBC activated with lipopolysaccharides (LPS), the anti-inflammatory properties of the studied extracts in comparison to piroxicam were investigated (Fig. 8B). It was inferred that Curcuma longa, Zingiber officinale, Salvia officinalis and Glycyrrhiza glabra extracts (0.19 ± 0.015, 1.57 ± 0.1, 2.31 ± 0.2 and 6.08 ± 0.5 μg/mL, respectively), demonstrated comparable effective anti-inflammatory concentrations (EAICs) www.nature.com/scientificreports/  -16   MAPK1, HMOX1, MMP2, CCND1, IL4, RAF1, MAPK3, EGF,  CDK2, TP53, PPARG, PTGER1, FOS, CYCS, CASP3, NQO1,  NOS2, CTNNB1, MTOR, PTGS2, JUN, PTEN, MMP9, PIM1,  MAPK8, NFE2L2, GSTP1, ESR1, AKT1, VEGFA, MYC   hsa05418  Fluid shear stress and atherosclerosis  18  1.21E-14   HMOX1, MMP2, MAPK14, ICAM1, TP53, VCAM1, NOS3,  FOS, NQO1, CTNNB1, JUN, MMP9, SRC, MAPK8, NFE2L2,  GSTP1, AKT1, VEGFA   hsa04926  Relaxin signaling pathway  15  1.08E-11  MAPK1, MMP2, MAPK14, RAF1, MAPK3, NOS3, FOS www.nature.com/scientificreports/ with piroxicam (1.21 ± 0.1 μg/mL), indicating the potential anti-inflammatory activities of these extracts. Realtime polymerase chain reaction (PCR) was used to assess the gene expression of four pro-inflammatory mediators (TNF-α, IL-1β, INF-γ and IL-6) in both untreated and lipopolysaccharide (LPS)-treated WBCs in order to ascertain the genetic basis of the anti-inflammatory effect (Fig. 8C-F). LPS increased the expression of the IL-1β gene by 6.3 ± 0.1-folds. When Allium sativum, Glycyrrhiza glabra, Oleu europea, Curcuma longa, Salvia officinalis L and Zingiber officinale were applied to the WBCs, this upregulation was eradicated by 0.5 ± 0.02, 0.3 ± 0.008, 3.6 ± 0.01, 0.4 ± 0.002, 1.3 ± 0.003 and 0.9 ± 0.003-fold, respectively to a degree that was comparable to that employed by piroxicam which exerted 0.4 ± 0.02-fold reduction in gene expression. In terms of INF-γ, LPS increased this gene expression by 7.3 ± 0.23-folds. This upregulation was reduced by 1.2 ± 0.05, 1.2 ± 0.03, 1.2 ± 0.009, 1.3 ± 0.01, and 1.5 ± 0.008-fold, respectively, after the WBCs were treated with Allium sativum, Salvia officinalis L, Glycyrrhiza glabra, Zingiber officinale and Curcuma longa extracts, respectively. This reduction was analogous to that caused by piroxicam (1.47 ± 0.03-fold). Meanwhile, LPS induced the TNF-α expression by 2.66 ± 0.01-folds that was diminished by the Glycyrrhiza glabra and Curcuma longa extracts to 0.016 ± 0.001 and 0.02 ± 0.03-folds, respectively. This level was nearly close to what piroxicam had achieved (0.03 ± 0.0003-folds). However, the expression of IL-6 was increased by LPS by 5.4 ± 0.2-folds, and this was reduced to 0.9 ± 0.02 and 1.3 ± 0.01-folds, respectively, by Allium sativum and Curcuma longa extracts. This amount was very comparable to what piroxicam had accomplished (1.1 ± 0.1-folds). It was noticed that Curcuma longa extract showed no statistically significant differences between it and the known synthetic anti-inflammatory drug, piroxicam. It's interesting to note that both the tested extracts and piroxicam considerably reduced the elevation of the gene expression of IL-1β, INF-γ, TNF-α and IL-6 to a similar degree (error bars were illustrated in Fig. 8 and all experiments p values were less than 0.01). Given the obvious suppression of the increased levels of TNF-α, IL-1β, INF-γ, IL-6 expression, it can be inferred that these studied extracts can serve as a possible natural antiinflammatory product. These findings were in line with network pharmacology studies, which revealed the tested extracts anti-inflammatory effect to have several targets and pathways.
Molecular docking studies of hit compounds in the active sites of the highest enriched target genes. The Glide module of the Schrodinger suite software was utilized to calculate the docking XP G scores of the top hit compounds apigenin, luteolin, diallyl trisulfide, silibinin, and allicin, in addition to piroxicam reference, against the active sites of the most enriched immunomodulatory target genes AKT1, CASP3, PTGS2, NOS3 and TP53. From Table 5, it can be noticed that silibinin had the smallest XP G score against RAC-alpha serine/threonine-protein kinase (AKT1) followed by caspase-3 and finally tumor suppressor protein TP53 with binding energies of − 11.725, − 8.855 and − 7.13 kcal/mol, respectively. Whereas luteolin and apigenin exhibited the highest stabilized interactions with RAC-alpha serine/threonine-protein kinase with binding energies of − 10.012 and − 9.013 kcal/mol, respectively followed by prostaglandin-endoperoxide synthase 2 with binding energies of − 7.534 and − 6.811 kcal/mol, respectively and then tumor suppressor protein TP53 with binding energies of − 6.732 and − 6.685 kcal/mol, respectively.
While the majority of interactions with the tumor suppressor protein TP53 (PDB ID 2VUK), in which apigenin was particularly integrated, were hydrophobic. Pro151, Pro223, Pro222, Cys220, Cys229, Trp146, Val147 and Leu145 were all found to interact hydrophobically with apigenin, and the 5, 4'-hydroxyl groups were particularly integrated in two hydrogen bonds with Cys220 and Asp228. This may have a significant impact on apigenin increased stability in the binding pocket of the target protein. Polar interactions with Thr150, Thr155, Thr230, and Ser149 were also present, in addition to charged negative interactions with Glu221, Asp228 and Asp148 (Fig. 9C). All of these investigations could successfully confirm the outcomes indicated above, suggesting the feasibility of this integrated strategy and serving as a guide to systematically uncover the therapeutic mechanisms of these naturally occurring plant components for the treatment of inflammation.

Molecular dynamics simulations stability test for hit protein-ligand complexes.
Two hundred (200) ns duration molecular dynamics (MD) simulation was performed to analyze the in silico stability of TP53 & apigenin, AKT1 & silibinin, and PTGS2 & luteolin protein-ligand complexes obtained from the molecular docking study 70 . The root mean square deviation (RMSD) measurements were performed for the shifts of ligands with respect to residues at the binding sites, and root mean square fluctuation (RMSF) measurements were made from the MD trajectory to explain fluctuations in protein structure 71 . As given in Fig. 10A Schematic 2D diagrams of protein-ligand interactions at 200 ns are given in Fig. 11 to demonstrate protein-ligand stabilities 72 . As shown in Fig. 11A, apigenin maintains its stability by forming an H bond with Leu145 at the active site of TP53, a negative charge with Asp228, and hydrophobic interactions with Cys220. As given in Fig. 11B, silibinin remained stable in the AKT1 active pocket, as in the molecular docking pose, giving Pi-Pi stacking interaction with Trp80, as well as an H bond interaction with Thr272 with Pi-Pi stacking and Cys296.  Table 3. Literature survey summary on the top scoring immunomodulatory plant constituents used as antiinflammatory agents. www.nature.com/scientificreports/ As given in Fig. 11C, luteolin gave Pi-Pi stacking interactions mainly with Thr355 at the PTGS2 active site and maintained its stability. Finally, animation videos from MD trajectory were created and provided in supplementary materials (Videos S1-S3) to monitor protein-ligand interactions full-time and better explain ligand stability. Analysis of MD simulations revealed that all three protein-ligand complexes are quite stable.

Experimental
The study design is summarized as shown in Fig. 12.
Assembling of an in-house database. Depending on our previous study of the chemical makeup of 32 selected immunomodulatory plants, 2154 phytochemical constituents were obtained for compilation of an in-house database as illustrated in Table S3. These plants were selected according to their traditionally known action to treat immune related disorders or based on prior in vivo, in vitro or clinical trials from literature review as mentioned in Khairy et al. 73 and summarized briefly in Table 6. Different resources including Dictionary of natural products (http:// dnp. chemn etbase. com/ faces/ chemi cal/ Chemi calSe arch. xhtml; jsess ionid= 0A16B F5251 5E734 B15A9 6DCBA E7788 B9), National Centre for Biotechnology Information's PubChem database (https:// pubch em. ncbi. nlm. nih. gov/), ChEMBL (https:// www. ebi. ac. uk/ chembl/) and Chemdraw software (CambridgeSoft Corporation, Cambridge, USA) were used to provide 2D structures of these compounds in (.sdf) format. Schrodinger software (2017A) was utilized to transform these 2D structures to SMILES format.   www.nature.com/scientificreports/ ADME and drug-likeness filtration. Using the Qikprop software (Schrodinger suite 2017A), plant constituents obtained from the database were filtered out through the application of ADME system and Lipinski's rule of five 28 . In order to determine if the examined compounds would be good therapeutic candidates, a variety of physio-chemical attributes were estimated. In this consider, compounds with an anticipated oral bioavailability (OB) ≤ 30 were not included in this investigation. In expansion, compounds that met fewer than three of Lipinski's five requirements were also disregarded.
Network pharmacology-based analysis. Target genes related to the filtered constituents. The target genes associated with the chosen chemical constituents were identified utilizing STITCH DB (http:// stitch. embl. de/, ver. 5.0) using the 'Homo sapiens' species setting 25 . Genes information like genes ID, names, organism and www.nature.com/scientificreports/   www.nature.com/scientificreports/ function were obtained from UniProt (http:// www. unipr ot. org/) 25,26 . Only the 'Homo sapiens' proteins connected to immunosuppressive disorders were kept. Then, using STRING database (https:// string-db. org) 36 , the protein-protein interaction network (PPI network) was created.  www.nature.com/scientificreports/ works including plant-constituent, gene-pathway, constituent-target gene and constituent-gene-pathway networks-were built using Cytoscape 3.8.2 (http:// www. cytos cape. org/). The nodes in these networks stand for genes, constituents and pathways, while the edges denote the interactions existing between them. The Cytoscape network analyzer plug-in was used to compute the network parameters. The importance of nodes in each built network was evaluated utilizing Cytoscape combined score of interactions.

Plant name Family Pharmacological action References
Allium sativum Amaryllidaceae Prevents the pro-inflammatory cytokines IL-6 and monocyte chemoattractant protein-1 (MCP-1) 74 Artemisia annua Asteraceae Suppresses the delayed-type hypersensitivity reaction-Had inhibitory effect on calmodulin 75,76 Calendula officinalis Asteraceae Possess anti-inflammatory activities-inhibits the mitogen-induced lymphocyte proliferation 77 Camellia sinensis Theaceae Possess anti-inflammatory properties 74 Cannabis sativa Cannabaceae Inhibits the proliferation of lymphocytes 78 Cichorium intybus Asteraceae Inhibits the lymphocyte proliferation assay in the presence of PHA (phytohemagglutinin) 77 Citrullus colocynthis Cucurbitaceae Reduces the cell proliferation induced by concanavalin A (con-A) 79 Citrus aurantium Rutaceae Active against concanavalin A and LPS (lipopolysaccharide) induced proliferation in both thymocytes and splenocytes, respectively 80

Propolis
Bee product High doses inhibit lymphocyte proliferation and possesses anti-inflammatory properties 96 Punica granatum Lythraceae Inhibits the activation of the nuclear factor of activated T cells, decreases CD3 + T cell infiltration of the inflamed tissue 75,97 Salvia officinalis L Lamiaceae Decreases the inflammation induced by LPS and reduces the blood levels of TNF-α and NF-κB 98

Sesamum indicum
Pedaliaceae Decreases white blood cell count, erythrocyte sedimentation rate, IL-6 and TNF-α 99 Silybum marianum Asteraceae Silymarin reduces Th1-linked cytokines (IL-2, IFN-γ, TNF-α) in addition to reduction of MAPKs' activities (ERK1/2 and P38) 100 Tanacetum parthenium Asteraceae Reduces several pro-inflammatory enzymes and mediators including phosphodiesterase-3 and 4, 5-lipoxygenase, NO, PGE2, IL-4 and TNF-α Genes and Genomes (KEGG) pathways and the Database for Annotation, Visualization and Integrated Discovery (DAVID) ver. 6.8 (https:// david. ncifc rf. gov/) were searched to learn more about gene ontology and to identify the canonical pathways, cellular components, biological processes and molecular functions that were closely related to the target genes 36,105 . The only selected pathways had P-values ≤ 0.05. P ≤ 0.05 was considered to be statistically significant.
Evaluation of network pharmacology analysis. The network pharmacology evaluation is conducted from three aspects-reliability, standardization, and rationality. Taking in consideration that the network pharmacology evaluation process includes both general evaluation and scalability evaluation. The evaluation process was carried out according to the method described previously by Li 68 .
Preparation of the extracts test solutions. According  Experimental verification of top target proteins. RAC-alpha serine/threonine-protein kinase (Akt1) inhibitory activity using spectrophotometric assay. The Akt1 inhibitory activity of the top four plants resulted from network pharmacology analysis was measured using a colorimetric assay of released adenosine diphosphate (ADP) 106 . The ADP assay kit offers a straightforward method for detecting ADP in a wide range of samples, specifically those containing reducing agents that may interfere with oxidase-based assays. ADP concentration was evaluated using a coupled enzyme colorimetric assay (450 nm) that was proportional to the quantity of ADP present. The bioassays were performed as described by Al-Sha' er et al. 107 and mentioned in supplementary material. Saturosporine, an Akt inhibitor, was tested as a positive control, while the negative controls were prepared by adding the substrate after terminating the reaction.
Caspase-3 inhibitory activity assay. The inhibitory activity of the top four plants resulted from network pharmacology analysis on caspase 3 gene was quantified as illustrated on caspase 3 colorimetric assay kit (CASP-3-C) 108 supplied by Sigma, and described by Ros et al. 109 as in supplementary materials.

Assessment of the cytotoxicity and anti-inflammatory activity.
The MTT test was used to evaluate the cytotoxicity of the various extracts in comparison to piroxicam. Each extract's effective anti-inflammatory concentrations (EAICs) were calculated in a culture of lipopolysaccharides (LPS) activated human WBCs obtained from human volunteers with ethical approval from ethical committee (IRB NO: 00012098). Expression levels of IL-6, IL-1β, TNF-α and INF-γ were determined using real-time polymerase chain reaction (PCR). The means ± standard deviations of three distinct replicates were used to express the results. The method described by Darwish et al. 110,111 , as given in the supplementary data, contains the specific details about the procedures. Statistical analysis of the results was performed using One-way ANOVA test-Tukey post-test options implemented on GraphPad Prism 5 program.

Molecular docking studies.
Crystallographic structure of the most enriched target genes (AKT1, CASP3, PTGS2, NOS3 and TP53), identified through network pharmacology analysis, were obtained from the Protein Data Bank (PDB). Each protein crystal structure was chosen based on the highest resolution possible. Through the use of the protein preparation wizard (OPLS 3 force field) module run in the Schrodinger suit, the preparation and energy minimization of crystal structures of the target proteins were carried out. The protein optimization was followed by the assignment of hydrogen bonds and bond order. At pH 7, zero order bonds to metals and disulphide bonds were also constructed. Next, all water molecules that were farther than 5°A away from the active site were eliminated. The grid boxes were constructed using the residues implicated in interactions with the co-crystallized ligands. The compounds 2D structures were loaded as (.sdf) format into the Lig Prep 2.3 module (Lig Prep, version 2.3, 2015, Schrödinger, USA) to build the least energy 3D structure for each compound and search for alternate conformers. The ionisation states were adjusted to create all conceivable states at pH 7. Molecular docking analyses were performed utilizing the Glide docking program of Maestro molecular modelling package executing extra-precision (XP-Glide) module. The ligand-target interactions such as hydrogen bond, hydrophobic interactions, ion pair interactions together with the binding modes of the identified compounds were demonstrated using Maestro interface. www.nature.com/scientificreports/

Molecular dynamics simulations.
To test the stability of protein-ligand complexes derived from molecular docking studies with Glide, MD simulations were run with Gromacs v2020.1 112 . The CHARMM-GUI server was used to prepare the necessary input files 113 . Using AMBER99SB force fields, proteins and compounds' topology files were created 114 . The TIP3 water model was used to dissolve the protein-ligand complexes, using a rectangular box type 10 Å away from the complexes. Then, 0.15 KCl salt was added to neutralize the solution. The created system was minimized to 5000 nsteps with the steep integrator. Nose-Hoover and Parrinello-Rahman algorithms were used to equilibrate it using NVT/NPT ensemble steps of 0.3 ns duration at 300 K and 1 atm pressure. 200 ns MD simulation was run and 1000 frame was recorded. The root mean square deviation (RMSD) and the root mean square fluctuation (RMSF) were measured with gmx scripts 115,116 . RMSD and RMSF plots were created with QtGrace v2.0.6. The PyMol Molecular Graphics System version 2.5.2 software was used to make MD trajectory videos and protein-ligand binding poses visualization.

Conclusion
Given that medicinal plants are multi-component and the network pharmacology approach emphasizes the idea of "network target, multicomponent therapeutics, " in a holistic manner analogous to the complex matrices of medicinal herbals, this strategy is thought to be appropriate for understanding the mechanism of action of medicinal plants. In this study, network pharmacology-based analysis of 2154 phytochemical constituents obtained from 32 selected immunomodulatory plants was carried out. A total of 37 constituents and 32 targets were found to be involved in 40 immunomodulatory-associated pathways. Apigenin, luteolin, diallyl trisulfide, silibinin and allicin had the highest percentage of C-T interactions, where the target genes AKT1, CASP3, PTGS2, NOS3, TP53 and MMP9 were found to be the most enriched ones by possessing the highest combined scores, suggesting that they may be the key nodes in C-T network. KEGG analysis illustrated that pathways in cancer, fluid shear stress and atherosclerosis, relaxin signaling pathway, IL-17 signaling pathway and FoxO signaling pathway had the largest number of observed genes and the smallest false discovery rate. Additionally, a combined plant-constituent-target-pathway network demonstrated that Curcuma longa, Allium sativum, Oleu europea, Salvia officinalis L, Glycyrrhiza glabra and Silybum marianum had the highest number of P-C-T-P interactions which would suggest that these plants have more active substances that can act as immunomodulatory agents. However, the main drawback in the current network pharmacology analysis studies is the insufficient scientific verification as it needs extensive in-vitro, in-vivo and clinical investigations 68 . Furthermore, molecular docking analysis of the top hit compounds; apigenin, luteolin, diallyl trisulfide, silibinin, and allicin, against the active sites of the most enriched immunomodulatory target genes; AKT1, CASP3, PTGS2, NOS3 and TP53 was carried out. The results revealed that silibinin had the lowest binding energy with RAC-alpha serine/threonine-protein kinase (AKT1) followed by caspase-3 and finally tumor suppressor protein TP53, whereas luteolin and apigenin exhibited the highest stabilized interactions with RAC-alpha serine/ threonine-protein kinase followed by prostaglandin-endoperoxide synthase 2 and then tumor suppressor protein TP53. The extracts of the highest scoring plants retrieved from network analysis, were then subjected to in vitro anti-inflammatory and cytotoxicity testing exhibiting outcomes that are equivalent to those of piroxicam. These plants are suggested as prospective sources of immunomodulatory agents by this study, which also offers a thorough explanation of the mechanism underlying the extracts' supposed anti-inflammatory action. To support our findings, additional in vivo and clinical research are advised.

Data availability
All data generated or analyzed during this study are included in this article (and its supplementary information files).